

/***************************************************************************************/ 
/*                                     _         _____                                 */
/*                       ___ ___ ___ _| |___ ___|_   _|___ _ _                         */
/*                      |  _| -_|   | . | -_|  _| | | | . | | |                        */
/*                      |_| |___|_|_|___|___|_|   |_| |___|_  |                        */
/*                                                        |___|                        */
/***************************************************************************************/


#include <stdio.h>
#include <ImathVec.h>

#include <ImfHeader.h>
#include <ImfRgbaFile.h>
#include <ImfArray.h>
#include <ImfStringAttribute.h>
#include <ImfMatrixAttribute.h>

#include <math.h>
#include <stdlib.h>
#include <float.h>

#include <pthread.h>
#include <vector>
#include <queue>
#include <algorithm>

#include <boost/shared_ptr.hpp>
#include <boost/lexical_cast.hpp>


using namespace std;
using namespace Imath;



const float epsilon = 10.f * FLT_EPSILON;


/***************************************************************************************/
/* Camera                                                                              */
/***************************************************************************************/


class Camera
{
public:  

    Camera(float near=0.01f, float fovy=60.f) : m_near(near), m_fovy(fovy)
    { 
        m_viewTarget = V3f(0);
		m_viewPoint = V3f(1);
        m_localX = V3f(1,0,0);
		m_localY = V3f(0,0,1);
        _update();
    }
    
    void set(V3f from, V3f target, V3f up) 
    {
        m_viewTarget = target;
        m_viewPoint = from;
        m_localY = up; 
        _update();
    }
    
    V3f getLocalX() const { return m_localX; }
    V3f getLocalY() const { return m_localY; }
    
    float m_fovy;
    float m_near;
    V3f m_viewPoint;
    V3f m_viewDirection;
    V3f m_viewTarget;
    V3f m_localX;
    V3f m_localY;
    
protected:
    
    void _update()
    {
        m_viewDirection = m_viewTarget - m_viewPoint;
        m_viewDirection = m_viewDirection.normalize();
        m_localX = m_viewDirection % m_localY;
        m_localY = m_localX % m_viewDirection;
        m_localX = m_localX.normalize();
        m_localY = m_localY.normalize();
    }
};



/***************************************************************************************/
/* Light                                                                               */
/***************************************************************************************/

class Light
{   
public:
    Light() : m_color(1,1,1), m_position(1,1,1) { }
    V3f m_color;
    V3f m_position;
};

typedef boost::shared_ptr<Light> LightPtr;


/***************************************************************************************/
/* Scene                                                                               */
/***************************************************************************************/

class Scene
{		
public:
    
    Scene(Box3f boundingBox, vector<LightPtr> lights) : 
        
        m_frame(0), m_maxAADepth(2), m_maxColorSep(0.03f), m_numRenderThreads(16), m_boundingBox(boundingBox), m_lights(lights) 
    
    {}
    
    inline void setFrame(int frame) { m_frame = frame; }
    inline int getFrame() { return m_frame; }
    
    inline void setMaxAADepth(int depth) { m_maxAADepth = depth; }
    inline int getMaxAADepth() { return m_maxAADepth; }
    
    inline void setMaxColorSep(float maxColorSep) { m_maxColorSep = maxColorSep; }
    inline float getMaxColorSep() { return m_maxColorSep; }
    
    inline void setNumRenderThreads(int numRenderThreads) { m_numRenderThreads = numRenderThreads; }
    inline int getNumRenderThreads() { return m_numRenderThreads; }
    
    // This gets a start point, and a normalized direction heading away from the eye.
    // It must simply return the color of the light returning to the eye along the opposite direction.
    virtual V3f shadePrimary(V3f start, V3f direction) const=0;
    
protected:
    
    int m_frame;
    Box3f m_boundingBox;
    
    // parameters driving adaptive anti-aliasation
    int m_maxAADepth;
    float m_maxColorSep;
    
    int m_numRenderThreads;
    
    vector<LightPtr> m_lights;
    
};


/***************************************************************************************/
/* Image                                                                               */
/***************************************************************************************/

struct Image
{
    Image(int width, int height) : m_width(width), m_height(height) 
    { 
        m_array = new Imf::Array2D<Imf::Rgba>(height, width);
    }
    
    ~Image()
    {
        delete m_array;
    }
    
    int width() { return m_width; }
    int height() { return m_height; }
    
    void setPixel(int i, int j, V3f color)
    {
        Imf::Rgba &p = (*m_array)[m_height-1-j][i];
        p.r = color.x;
        p.g = color.y;        
        p.b = color.z;
        p.a = 1.0;
    }
    
    void writeToFile(const char fileName[])
    {
        const Imf::Rgba *pixels = &(*m_array)[0][0];
        Imf::RgbaOutputFile file (fileName, m_width, m_height, Imf::WRITE_RGBA);
        file.setFrameBuffer (pixels, 1, m_width);
        file.writePixels (m_height);
    }
    
    private:
        
        int m_width;
        int m_height;
        Imf::Array2D<Imf::Rgba>* m_array;

};


/***************************************************************************************/
/* Raycaster                                                                           */
/***************************************************************************************/

struct Ray
{
    V3f m_origin;
    V3f m_direction;
};

// Four raycasts at the corners of a square form a chunk. Multiple chunks are combined
// to render a single antialiased pixel.
struct Chunk
{
    Chunk() {}
    Chunk(float L, float cW, float cH, V3f TL, V3f TR, V3f BL, V3f BR) :
        m_L(L), m_cW(cW), m_cH(cH), m_TL(TL), m_TR(TR), m_BL(BL), m_BR(BR) {}
    
    V3f m_TL, m_TR, m_BL, m_BR;
    float m_cW, m_cH;
    float m_L;
    
    inline V3f getAverage() const
    {
        V3f C;
		for (int a=0; a<3; ++a) C[a] = m_BL[a] + m_BR[a] + m_TL[a] + m_TR[a];
		return C/4.0f;
    }
    
    inline float getSeparation(const V3f& average) const
    {
        V3f tl_diff = m_TL - average; 
		V3f tr_diff = m_TR - average; 
		V3f bl_diff = m_BL - average; 
		V3f br_diff = m_BR - average; 
        float tl = tl_diff ^ tl_diff; 
		float tr = tr_diff ^ tr_diff;
		float bl = bl_diff ^ bl_diff;
		float br = br_diff ^ br_diff;
		return max(max(tl, tr), max(bl, br));
    }
};

struct Task
{
    Task(int X, int Y, int W, int H, float invFW, float invFH, float frustum_width, float frustum_height, const Camera& camera, Scene& scene, Image& image) :
    
        m_X(X), m_Y(Y), m_W(W), m_H(H), m_invFW(invFW), m_invFH(invFH), m_frustum_width(frustum_width), m_frustum_height(frustum_height), 
        m_camera(camera), m_scene(scene), m_image(image)
    
    {
        // Storage of previous and current task line
		m_lastLine.resize(m_W+1);
		m_newLine.resize(m_W+1);
		m_line.resize(m_W+1);
    }
    
    ~Task()
    {
        m_lastLine.clear();
		m_newLine.clear();
		m_line.clear();
    }
    
    void run()
    {
        for (int j = m_X; j<= m_X + m_W; ++j)		
        {
            Ray primary = makePrimaryRay(float(j), float(m_Y));
            m_lastLine[j-m_X] = m_scene.shadePrimary(primary.m_origin, primary.m_direction);
        }
        
        // Do the rest of the task rectangle
		for (int i=m_Y-1; i>=m_Y-m_H; --i)
		{
            Ray primary = makePrimaryRay(float(m_X), float(i));
            
            // Each primary chunk uses raycasts from the previous pixel and line.
			// Provide the previous value for the initial pixel of each line.
            m_lastBR = m_scene.shadePrimary(primary.m_origin, primary.m_direction);
            for (int j = 0; j<m_W; ++j)			
			{
                // Each pixel has an associated "primary chunk" of 4 raycasts through its corners
				// (although only one new primary raycast needs to be done for each pixel)
				float cW = m_X + j + 0.5f;
				float cH = i + 0.5f;
                Chunk chunk = makePrimaryChunk(cW, cH, i, j);
                
                // Subdivide primary chunk into sub-chunks for adaptive anti-aliasing.
				V3f color = adaptiveAntialiasation(chunk, 0);
                //expose(color, exposure);
				//gammaCorrect(color);
                m_line[j] = color;
            }
            
            // Plot line to framebuffer
			for (int j = 0; j<m_W; ++j)			
			{
				const V3f& color = m_line[j];
                m_image.setPixel(m_X + j, i, color);
			}
            
            for (int i=0; i<m_W+1; ++i) m_lastLine[i] = m_newLine[i];
        }
    }
    
    inline Ray makePrimaryRay(float pW, float pH)
    {
        Ray ray;
        ray.m_origin = m_camera.m_viewPoint;
        float lW = 0.5f * m_frustum_width  * (2.0f * pW * m_invFW - 1.0f);
		float lH = 0.5f * m_frustum_height  * (2.0f * pH * m_invFH - 1.0f);
        V3f O = m_camera.m_near * m_camera.m_viewDirection;
        V3f X = m_camera.getLocalX();
        V3f Y = m_camera.getLocalY();
        V3f dir = O + lW*X + lH*Y;
        ray.m_direction = dir.normalize();
        return ray;
    }
    

    inline Chunk makePrimaryChunk(float cW, float cH, int i, int j)
    {
        Chunk chunk;
        chunk.m_cW = cW;
        chunk.m_cH = cH;
		chunk.m_L = 1.0f;
		
		//top left
		chunk.m_TL = m_lastLine[j];
		
		//top right
		chunk.m_TR = m_lastLine[j+1];
		
		//bottom left
		chunk.m_BL = m_lastBR;
		m_newLine[j] = chunk.m_BL;
        
        //bottom right (the only new primary ray cast per pixel)
		Ray primary = makePrimaryRay(cW + 0.5f, cH - 0.5f);
        chunk.m_BR = m_scene.shadePrimary(primary.m_origin, primary.m_direction);
		m_newLine[j+1] = chunk.m_BR;
		
		//keep BR color for next pixel BL
		m_lastBR = chunk.m_BR;
		return chunk;
    }
    
    
    inline V3f adaptiveAntialiasation(Chunk chunk, int depth)
    {
        V3f average = chunk.getAverage();
        if ( m_scene.getMaxAADepth() == depth ) return average;
        
        // compute maximum color separation from average
		float sep = chunk.getSeparation(average);
		if (sep < m_scene.getMaxColorSep() * m_scene.getMaxColorSep()) return average;
        
        // construct 4 sub-chunks and recurse on each
		depth++;
        V3f TOP = chunkCast(0, 1, chunk); 
		V3f BOT = chunkCast(0, -1, chunk); 
		V3f LEF = chunkCast(-1, 0, chunk); 
		V3f RIH = chunkCast(1, 0, chunk); 
		V3f CEN = chunkCast(0, 0, chunk);

        float hL = 0.5f  * chunk.m_L; 
		float qL = 0.25f * chunk.m_L;
		float cW_neg = chunk.m_cW - qL;
		float cW_pos = chunk.m_cW + qL;
		float cH_neg = chunk.m_cH - qL;
		float cH_pos = chunk.m_cH + qL;
		Chunk subTL(hL, cW_neg, cH_pos, chunk.m_TL, TOP, LEF, CEN);
		Chunk subTR(hL, cW_pos, cH_pos, TOP, chunk.m_TR, CEN, RIH);
		Chunk subBL(hL, cW_neg, cH_neg, LEF, CEN, chunk.m_BL, BOT);
		Chunk subBR(hL, cW_pos, cH_neg, CEN, RIH, BOT, chunk.m_BR);
        
        V3f colTL = adaptiveAntialiasation(subTL, depth);
		V3f colTR = adaptiveAntialiasation(subTR, depth);
		V3f colBL = adaptiveAntialiasation(subBL, depth);
		V3f colBR = adaptiveAntialiasation(subBR, depth);
        Chunk smooth(hL, chunk.m_cW, chunk.m_cH, colTL, colTR, colBL, colBR);
		return smooth.getAverage();
    }


    inline V3f chunkCast(int i, int j, const Chunk& chunk)
    {
        float pW = chunk.m_cW + i * 0.5f * chunk.m_L;
		float pH = chunk.m_cH + j * 0.5f * chunk.m_L;
		Ray primary = makePrimaryRay(pW, pH );
		return m_scene.shadePrimary(primary.m_origin, primary.m_direction);
    }
    
    // Coordinates of the top-left pixel of this task in the image
    int m_X;
    int m_Y; 
    
    // Pixel dimensions of task region (square except at edge) 
    int m_W;
    int m_H;
    
    // Inverse width and height of entire framebuffer
    float m_invFW;
    float m_invFH;
    
    // Camera frustum dimensions in world space lengths
    float m_frustum_width;
    float m_frustum_height;
    
    const Camera& m_camera;
    Scene& m_scene;
    Image& m_image;
    
    // Storage
    vector<V3f> m_lastLine;
    vector<V3f> m_newLine;
    vector<V3f> m_line;
    V3f m_lastBR;
    
};


// Rendering is broken into tasks run in parallel, each task working on a different section of the framebuffer.

struct TaskPool
{
    TaskPool() : m_renderAborted(false)
    {
        const pthread_mutexattr_t* attr = NULL;
        int err = pthread_mutex_init(&m_taskMutex, attr);
        if (err != 0)
        {
            fprintf(stderr, "Could not pthread_mutex_init, error %d. %s\n", err, strerror(err));
            exit(EXIT_FAILURE);
        }
    }
    
    std::deque<Task> m_taskQueue;
    pthread_mutex_t m_taskMutex;
    bool m_renderAborted;
};


struct Thread
{
    pthread_t m_threadID;
};


static void* _threadFunc(void* data)
{
    TaskPool* taskPool = reinterpret_cast<TaskPool*>(data);
    while (!taskPool->m_renderAborted)
    {
        // Enter critical section while trying to pull a new task from the queue
        pthread_mutex_lock(&taskPool->m_taskMutex);
        int nRemainingTasks = taskPool->m_taskQueue.size();
        
        // If the task queue is empty, the thread quits
        if ( nRemainingTasks==0 )
        {
            pthread_mutex_unlock(&taskPool->m_taskMutex);
            return NULL;
        }
        
        // Pull a task from the queue
        Task task = taskPool->m_taskQueue.back();
        taskPool->m_taskQueue.pop_back();
        pthread_mutex_unlock(&taskPool->m_taskMutex);
        
        // Do the rendering task
        task.run();
    }
    
    return NULL;
}


void render(Image& image, Scene& scene, const Camera& camera)
{
    TaskPool m_taskPool;
    
    int FW = image.width();
    int FH = image.height();
    float invFW = 1.f/float(FW);
    float invFH = 1.f/float(FH);
    float aspect = float(FW)/float(FH);
    
    // (FW, FH) = frustum dimensions
    float frustumHeight = 2.0f * camera.m_near * tan(camera.m_fovy * M_PI/360.f);
    float frustumWidth = frustumHeight * aspect;
    const int TX = 32;
    const int TY = 32;
    int numTaskX = ( FW % TX == 0) ?  FW / TX : ceil(float(FW)/float(TX));
    int numTaskY = ( FH % TY == 0) ?  FH / TY : ceil(float(FH)/float(TY));
    
    for (int i=0; i<numTaskX; ++i)
    {
        for (int j=0; j<numTaskY; ++j)
        {
            int X = i * TX;
            int Y = j * TY;
            int W = ( FW - X >= TX ) ? TX : FW - X;
            int H = ( FH - Y >= TY ) ? TY : FH - Y;
            Task task(X, FH-Y, W, H, invFW, invFH, frustumWidth, frustumHeight, camera, scene, image);                 
            m_taskPool.m_taskQueue.push_back(task);
        }
    }
    
    // A set number of threads pull tasks from the queue
    vector<Thread> threads;
    
    for (int nThread=0; nThread<scene.getNumRenderThreads(); ++nThread)
    {
        // Create the thread using POSIX routines.
        pthread_attr_t  attr;
        assert(!pthread_attr_init(&attr));
        
        Thread thread;
        void* data = (void*) &m_taskPool;
        int threadError = pthread_create(&thread.m_threadID, &attr, &_threadFunc, data);
        assert(!pthread_attr_destroy(&attr));
        if (threadError != 0)
        {
            fprintf(stderr, "Could not pthread_create thread %d. %d/%s\n", nThread, threadError, strerror(threadError));
            exit(EXIT_FAILURE);
        }
        threads.push_back(thread);
    }
    
    cout << "waiting for render tasks to finish ..." << endl;
    
    // Wait until all threads are done
    for (int nThread=0; nThread<scene.getNumRenderThreads(); ++nThread)
    {
        Thread& thread = threads[nThread];
        void* retVal;
        int threadError = pthread_join(thread.m_threadID, &retVal);
        if (threadError != 0)
        {
            fprintf(stderr, "Error joining thread %d. %d/%s\n", nThread, threadError, strerror(threadError));
            exit(EXIT_FAILURE);
        }
    }
    
    cout << "Render tasks complete." << endl;
}





/***************************************************************************************/
/*  Example calling code                                                               */
/***************************************************************************************/


/******************************************************************************************************************/
// Perlin noise code

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define B 0x100
#define BM 0xff
#define Nperlin 0x1000
#define NM 0xfff
#define s_curve(t) ( t * t * (3. - 2. * t) )
#define lerp(t, a, b) ( a + t * (b - a) )
#define setup(i,b0,b1,r0,r1)\
t = vec[i] + Nperlin;\
b0 = ((int)t) & BM;\
b1 = (b0+1) & BM;\
r0 = t - (int)t;\
r1 = r0 - 1.;
#define at2(rx,ry) ( rx * q[0] + ry * q[1] )
#define at3(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] )

static int p[B + B + 2];
static double g3[B + B + 2][3];
static double g2[B + B + 2][2];
static double g1[B + B + 2];
static int start = 1;

namespace PerlinNoise
{
    void init(void);
    double noise1(double);
    double noise2(double *);
    double noise3(double *);
    void normalize3(double *);
    void normalize2(double *);
    
    double PerlinNoise1D(double,double,double,int);
    double PerlinNoise2D(double,double,double,double,int);
    double PerlinNoise3D(double,double,double,double,double,int);
    
    double noise1(double arg)
    {
        int bx0, bx1;
        double rx0, rx1, sx, t, u, v, vec[1];
        
        vec[0] = arg;
        if (start) {
            start = 0;
            init();
        }
        
        setup(0,bx0,bx1,rx0,rx1);
        
        sx = s_curve(rx0);
        u = rx0 * g1[ p[ bx0 ] ];
        v = rx1 * g1[ p[ bx1 ] ];
        
        return(lerp(sx, u, v));
    }
    
    double noise2(double vec[2])
    {
        int bx0, bx1, by0, by1, b00, b10, b01, b11;
        double rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
        int i, j;
        
        if (start) {
            start = 0;
            init();
        }
        
        setup(0, bx0,bx1, rx0,rx1);
        setup(1, by0,by1, ry0,ry1);
        
        i = p[ bx0 ];
        j = p[ bx1 ];
        
        b00 = p[ i + by0 ];
        b10 = p[ j + by0 ];
        b01 = p[ i + by1 ];
        b11 = p[ j + by1 ];
        
        sx = s_curve(rx0);
        sy = s_curve(ry0);
        
        q = g2[ b00 ] ; u = at2(rx0,ry0);
        q = g2[ b10 ] ; v = at2(rx1,ry0);
        a = lerp(sx, u, v);
        
        q = g2[ b01 ] ; u = at2(rx0,ry1);
        q = g2[ b11 ] ; v = at2(rx1,ry1);
        b = lerp(sx, u, v);
        
        return lerp(sy, a, b);
    }
    
    double noise3(double vec[3])
    {
        int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
        double rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
        int i, j;
        
        if (start) {
            start = 0;
            init();
        }
        
        setup(0, bx0,bx1, rx0,rx1);
        setup(1, by0,by1, ry0,ry1);
        setup(2, bz0,bz1, rz0,rz1);
        
        i = p[ bx0 ];
        j = p[ bx1 ];
        
        b00 = p[ i + by0 ];
        b10 = p[ j + by0 ];
        b01 = p[ i + by1 ];
        b11 = p[ j + by1 ];
        
        t  = s_curve(rx0);
        sy = s_curve(ry0);
        sz = s_curve(rz0);
        
        q = g3[ b00 + bz0 ] ; u = at3(rx0,ry0,rz0);
        q = g3[ b10 + bz0 ] ; v = at3(rx1,ry0,rz0);
        a = lerp(t, u, v);
        
        q = g3[ b01 + bz0 ] ; u = at3(rx0,ry1,rz0);
        q = g3[ b11 + bz0 ] ; v = at3(rx1,ry1,rz0);
        b = lerp(t, u, v);
        
        c = lerp(sy, a, b);
        
        q = g3[ b00 + bz1 ] ; u = at3(rx0,ry0,rz1);
        q = g3[ b10 + bz1 ] ; v = at3(rx1,ry0,rz1);
        a = lerp(t, u, v);
        
        q = g3[ b01 + bz1 ] ; u = at3(rx0,ry1,rz1);
        q = g3[ b11 + bz1 ] ; v = at3(rx1,ry1,rz1);
        b = lerp(t, u, v);
        
        d = lerp(sy, a, b);
        
        return lerp(sz, c, d);
    }
    
    void normalize2(double v[2])
    {
        double s;
        
        s = sqrt(v[0] * v[0] + v[1] * v[1]);
        v[0] = v[0] / s;
        v[1] = v[1] / s;
    }
    
    void normalize3(double v[3])
    {
        double s;
        
        s = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
        v[0] = v[0] / s;
        v[1] = v[1] / s;
        v[2] = v[2] / s;
    }
    
    void init(void)
    {
        int i, j, k;
        
        for (i = 0 ; i < B ; i++) {
            p[i] = i;
            g1[i] = (double)((random() % (B + B)) - B) / B;
            
            for (j = 0 ; j < 2 ; j++)
                g2[i][j] = (double)((random() % (B + B)) - B) / B;
            normalize2(g2[i]);
            
            for (j = 0 ; j < 3 ; j++)
                g3[i][j] = (double)((random() % (B + B)) - B) / B;
            normalize3(g3[i]);
        }
        
        while (--i) {
            k = p[i];
            p[i] = p[j = random() % B];
            p[j] = k;
        }
        
        for (i = 0 ; i < B + 2 ; i++) {
            p[B + i] = p[i];
            g1[B + i] = g1[i];
            for (j = 0 ; j < 2 ; j++)
                g2[B + i][j] = g2[i][j];
            for (j = 0 ; j < 3 ; j++)
                g3[B + i][j] = g3[i][j];
        }
    }
    
    double PerlinNoise1D(double x,double alpha,double beta,int n)
    {
        int i;
        double val,sum = 0;
        double p,scale = 1;
        
        p = x;
        for (i=0;i<n;i++) {
            val = noise1(p);
            sum += val / scale;
            scale *= alpha;
            p *= beta;
        }
        return(sum);
    }
    
    double PerlinNoise2D(double x,double y,double alpha,double beta,int n)
    {
        int i;
        double val,sum = 0;
        double p[2],scale = 1;
        
        p[0] = x;
        p[1] = y;
        for (i=0;i<n;i++) {
            val = noise2(p);
            sum += val / scale;
            scale *= alpha;
            p[0] *= beta;
            p[1] *= beta;
        }
        return(sum);
    }
    
    double PerlinNoise3D(double x,double y,double z,double alpha,double beta,int n)
    {
        int i;
        double val,sum = 0;
        double p[3],scale = 1;
        p[0] = x;
        p[1] = y;
        p[2] = z;
        for (i=0;i<n;i++) {
            val = noise3(p);
            sum += val / scale;
            scale *= alpha;
            p[0] *= beta;
            p[1] *= beta;
            p[2] *= beta;
        }
        return(sum);
    }
}

#undef B
#undef BM
#undef Nperlin
#undef NM
#undef s_curve
#undef lerp
#undef setup
#undef at2
#undef at3

/******************************************************************************************************************/


 template<typename T>
 inline Vec3<T> vabs(Vec3<T> v)
 {
 return Vec3<T>(fabs(v.x), fabs(v.y), fabs(v.z));
 }
 

 template<typename T>
 inline Vec3<T> vmod(Vec3<T> v, T m)
 {
 return Vec3<T>(fmod(v.x, m), fmod(v.y, m), fmod(v.z, m));
 }
 

inline V3d vabs(V3d v)
{
    return V3d(fabs(v.x), fabs(v.y), fabs(v.z));
}

inline V3d vmod(V3d v, double m)
{
    return V3d(fmod(v.x, m), fmod(v.y, m), fmod(v.z, m));
}



/* copied from the code shown in the video here: http://www.subblue.com/blog/2011/3/5/fractal_lab */
inline double mengerSponge(V3d w, double scale=5.0, int maxIterations=4)
{
    V3d halfSpongeScale = V3d(0.5) * scale;
    
    w = (w*0.5 + V3d(0.5)) * scale;
    V3d v = vabs(w - halfSpongeScale) - halfSpongeScale;
    
    double dl = max(v.x, max(v.y, v.z));
    double d = dl;
    double p = 1.0;
    for (int i=0; i<maxIterations; ++i)
    {
        V3d a = vmod(3.0 * w * p, 3.0);
        p *= 3.0;
        v = V3d(0.5) - vabs(a - V3d(1.5));
        
        // distance inside the 3 axis-aligned square tubes
        dl = min(max(v.x, v.z), min(max(v.x, v.y), max(v.y, v.z))) / p;
        
        // intersection
        d = max(d, dl);
    }
    
    return d*2.0/scale;
}


class DistanceSurface: public Scene
{
    public:
    
        DistanceSurface(Box3f boundingBox, vector<LightPtr> lights) : 
    
            Scene(boundingBox, lights) 
    
        {
            m_maxIterations = 128;
            m_surfaceTol = 1.0e-6 * (m_boundingBox.max - m_boundingBox.min).length();
            m_epsilonDisplace = 1.5f*m_surfaceTol;
            m_invEpsilonDistance = 0.5 / m_epsilonDisplace;
            m_skyColor = 0.001*V3d(0.01,0.01,0.05);//0.01 * V3d(0.5,0.5,0.6);
        }
        
    
        int m_maxIterations;
        double m_surfaceTol; 
        double m_epsilonDisplace; 
        double m_invEpsilonDistance;
    
        V3d m_skyColor;
    
        
        inline float subtract(float d1, float d2) const
        {
            return fmax(-d1, d2);
        }
    

        inline double sphere(V3f P, double radius) const
        {
            return P.length() - radius;
        }
    
        inline double roundBox(V3f P, V3f edgeLengths, double radius) const
        {            
            V3d M(fmax(fabs(P.x) - edgeLengths.x, 0), 
                  fmax(fabs(P.y) - edgeLengths.y, 0), 
                  fmax(fabs(P.z) - edgeLengths.z, 0));
            return M.length() - radius;
        }
    
        inline double distanceToSurface(V3f P) const
        {
            return mengerSponge(P);
        }
    
        virtual V3f shadePrimary(V3f start, V3f direction) const
        {
            V3d P = start;
            bool hit = false;
            int nIterations = 0;
            
            V3d X(1,0,0);
            V3d Y(0,1,0);
            V3d Z(0,0,1);
            V3d normal;
            
            while (nIterations < m_maxIterations && P.length()<1000.0)
            {
                double unboundingRadius = distanceToSurface(P);
                if (unboundingRadius < m_surfaceTol)
                {
                    hit = true;
                    
                    // Compute normal by central differences
                    V3d pd[6];
                    pd[0] = P + m_epsilonDisplace*X;
                    pd[1] = P - m_epsilonDisplace*X;
                    pd[2] = P + m_epsilonDisplace*Y;
                    pd[3] = P - m_epsilonDisplace*Y;
                    pd[4] = P + m_epsilonDisplace*Z;
                    pd[5] = P - m_epsilonDisplace*Z;
                    
                    double gradX = ( distanceToSurface(pd[0]) - distanceToSurface(pd[1]) ) * m_invEpsilonDistance;
                    double gradY = ( distanceToSurface(pd[2]) - distanceToSurface(pd[3]) ) * m_invEpsilonDistance;
                    double gradZ = ( distanceToSurface(pd[4]) - distanceToSurface(pd[5]) ) * m_invEpsilonDistance;
                    normal.setValue(gradX, gradY, gradZ);
                    normal.normalize();															  
                    break;
                }
                
                // March along ray
                P += 0.99 * unboundingRadius * V3d(direction);
                nIterations++;
            }
            
            if (!hit)
            {
                return m_skyColor;
            }
            
            V3f color(0);
            for (int l=0; l<m_lights.size(); ++l)
            {
                LightPtr L = m_lights[l];
                V3f lightDir = (L->m_position - P).normalize();  
                
                float diffuse = fabs(normal^lightDir);
                float lambert = fmax(0.0, diffuse);
                
                V3f H = lightDir - 2.0f*diffuse*V3f(normal);
                H = H.normalize();
                float specPower = 60.f;
                float phong = powf(fabs(direction^H), specPower);
                
                color += L->m_color * (10.0*phong) + 0.1 * m_skyColor*lambert;
                //color += L->m_color * (phong);

            }
            return color;
        }
    
};


int main(int argc, const char * argv[])
{
    Camera camera(0.01, 60.f);
    V3f from(2.5, 2, 1.5), target(0, 0, 0), up(0, 0, 1);
    camera.set(from, target, up);
  
    std::vector<LightPtr> lights;
    {
        LightPtr l = LightPtr(new Light());
        l->m_color = V3f(0.7,0.6,0.3);
        l->m_position = V3f(500, 500, 0);
        lights.push_back(l);
    }
    {
        LightPtr l = LightPtr(new Light());
        l->m_color = V3f(0.1,0.1,0.1);
        l->m_position = V3f(-500, -500, 0);
        lights.push_back(l);
    }
    {
        LightPtr l = LightPtr(new Light());
        l->m_color = V3f(0.1,0.4,0.2);
        l->m_position = V3f(0, 500, 500);
        lights.push_back(l);
    }
    {
        LightPtr l = LightPtr(new Light());
        l->m_color = V3f(0.8,0.4,0.8);
        l->m_position = V3f(0, -500, -500);
        lights.push_back(l);
    }
    
    
    Box3f boundingBox;
    boundingBox.min = V3f(-1.0);
    boundingBox.max = V3f(1.0);
    
    DistanceSurface scene(boundingBox, lights);
    scene.setFrame(0);
    
    scene.setMaxAADepth(2);
    scene.setMaxColorSep(0.005);
    
    int width = 1980; 
    int height = 1080;
    Image image(width, height);
    
    for (int n=0; n<1; ++n)
    {
        scene.setFrame(n);
        camera.set(from, target, up);
        
        cout << "\n----------------------------------------" << endl;
        cout << "Rendering frame " << n << " ..." << endl;
        render(image, scene, camera);
        cout << "Done." << endl;
        
        string filename = string("/Users/jamports/projects/renderToy/RenderToy/testImages/roundbox");
        filename += string(".") + boost::lexical_cast<string>(n) + string(".exr");
        cout << "Writing image " << filename << " ..." << endl;
        image.writeToFile(filename.c_str());
        
        cout << "Done." << endl;
    }
    
    // then:  djv_convert movie0-100.exr output.mov
    
        
    return 0;
}







